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We present a new tool, GPA, that can generate key performance measures for very large systems. 
Based on solving systems of ordinary differential equations (ODEs), this method of performance 
analysis is far more scalable than stochastic simulation. The GPA tool is the first to produce higher 
moment analysis from differential equation approximation, which is essential, in many cases, to 
obtain an accurate performance prediction. We identify so-called switch points as the source of 
error in the ODE approximation. We investigate the switch point behaviour in several large models 
and observe that as the scale of the model is increased, in general the ODE performance prediction 
improves in accuracy. In the case of the variance measure, we are able to justify theoretically that in 
the limit of model scale, the ODE approximation can be expected to tend to the actual variance of the 
model. 

1 Introduction 

Quantitative analysis of systems by means of differential equation (ODE) approximation lUl IH or fluid 
techniques ||3l IH produce transient performance measures for massive state-space process models of 
lO'*"' states and beyond. Previous explicit state-space performance techniques which analysed the un- 
derlying continuous-time Markov chain directly (for example, EllH) were limited to 10^' states in the 
very best cases and then only for the simplest steady-state style of analysis. 

In physical and biological processes, deterministic approximation of system evolution via systems of dif- 
ferential equations have existed for some time ||71|8l|9l. However, differential equation-based techniques 
have only recently been brought to bear on process models of computer systems. A major difference 
between the two lies in the model of interaction assumed in the two distinct fields. The model of syn- 
chronisation used in computer and communication systems differs from that typically used in physical 
and biological processes, where for example mass-action dynamics govern the system evolutionj^ This 
difference significantly changes the nature of the differential equation approximation of computer sys- 
tems, thus results from mass-action-based systems cannot be translated to systems based on, for example, 
process algebras, queueing networks or stochastic Petri nets. 

In particular, there is not as yet a detailed understanding of how good the differential equation approx- 
imation is to the underlying discrete Markov process as generated from process models of computer 
systems. For instance [3] produced a first order approximation to large scale Markov models, but there is 
no discussion of accuracy of the approximation or higher moment generation. The issue focuses around 
so-called switch points in a model. In the case of PEPA, and other computer performance modelling 
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formalisms, these were identified as the source of the error in the differential equation approximations 
in im. We use our new GPA tool, to explore these switch points and the error associated with them. 

We aim to show that in many cases the existence and location of an approximation error in the fluid 
model can be predicted. Our tool will be able to warn the modeller about the presence of error for 
certain initial parameter regimes or in particular time-intervals for a given model execution. We examine 
not only first moment predictions for several simple performance models but also variances and higher 
moment solutions of the ODE approximation. In all cases it is essential to know when the analysis is 
accurate. Higher moments, in particular, can be used to create passage-time analysis bounds [10.1 and 
accurate moment approximations will make for precise bounds. 



1.1 Motivating Example 

We will first look at a simple motivating example. Consider the ubiquitous situation of m identical 
processors running in parallel, each in need of one of n identical resources. Each processor repeatedly 
has to acquire an available resource, after which it is ready to perform a required task and return to 
the initial state. Each resource has to be reset after it is acquired. The actions the system takes (e.g. a 
processor acquiring a resource or a resource resetting) are stochastic in nature and do not have a fixed 



duration. Formally, this model is defined in Section 1.3 



We are interested in how the system evolves over time: that is, in the counts of the four different states 
within the system (available and acquired resources and waiting and ready processors) at time, t. One 
way to do this is to simulate the system repeatedly and take the mean of counts of each state at each 
point of time. An efficient alternative is to deterministically approximate the evolution of the model 
by a system of coupled ordinary differential equations using techniques from, for example |[T]|2j|3]lll. 
Figure [T] shows a plot for each of these methods for our example, for m = 50 and n = 20. 



Waiting 
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(a) ODE approximation: mean 



Time, t 

(b) Simulation: mean 



Figure 1: Comparison of the simulation averages and numerical solution to the corresponding approximat- 
ing ODEs for the processor/resource model. 



It can be useful to know not just the average counts of states at each time but also the variability of these 
counts. For the stochastic simulation, this can be achieved by calculating the variance instead of just 
the mean. For the ODE approximation, it has been shown in how to adapt the existing techniques to 
produce systems of ODEs approximating higher moments of state counts. Figure[2]compares a numerical 
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solution to these ODEs with the variance obtained from the stochastic simulation. 



Waiting processor - 
Ready processor 
Avaiiable resource - 
Acquired resource 



Waiting p 
Ready processor 
Available r( 
Acquired rr 




(a) ODE approximation: variance 



(b) Simulation: variance 



Figure 2: Comparison of tlie variance for tlie processor/resource model obtained from tlie stocliastic simu- 
lation with tlieir ODE approximation. 



It is apparent that the ODE approximation is not equally accurate for the whole observed time. Figure [T] 
and even more so Figure [2] shows that the error of the approximation starts just before the time t = 0.2. 

The main focus of our work is to investigate, with the help of an efficient software tool, why this error 
appears and to provide practical means of detecting it without running the computationally expensive 
simulation. We will explain, using results from [1], that it occurs due to models passing through so- 
called switch points, when the total rate of cooperating actions between groups of components becomes 
equal. 

We present a tool that will visualise the dis- 
tance that an evolving model is from a switch 
point. Figure [3] shows this for the processor 
and resource model and confirms that the error 
in the variance approximation coincides with 
the switch point in total rate between proces- 
sors and resources just before t = 0.2. 

In this work, we will consider models described 
in an extension of the stochastic process al- 



gebra PEPA. Section 1.3 provides a brief in 



troduction to Grouped PEPA. In Section 1.3.1 



we will overview the method of deriving ODE 
approximations to higher moments of PEPA 
models. 

It was shown in flO) that the error of this ap- 
proximation for the mean decreases as we in- 
crease the initial populations (e.g. m and n in 




Figure 3: Switch point plot showing the difference be- 
tween total rate of processors and resources component 
groups. The switch point occurs when the plot goes 
through zero. 



this example). Section 1.5 gives a theoretical justification that a similar convergence can be expected to 
hold in the case of the variance ODE approximation. 
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It is important to note that although we are studying switch points in the context of PEPA models, switch 
points also occur in other performance modelling settings. Specifically, they occur whenever there is a 
contention for a finite pool of resources which may be saturated by demand. The most obvious example 
of this is an M/M/n queue which has n servers, and a service rate of jJ. at each server; the total service 
rate when k jobs are in the system is given by min(?i/x jkjj.), where k = n defines the switch point of this 
queue. 

Section |2] introduces a tool that for the first time implements higher order moment approximation for a 
stochastic process algebra. This tool is used for our investigations in Section |3j where we present cases 



studies of different switch point behaviour and also demonstrate instances of the results from Section 1 .5 



1.2 Related Tools 

There are many tools that support analysis of very large state spaces in performance modelling. Two 
such popular tools which have good support for explicit state-space analysis are Mobius and PRISM. 

The Mobius fTTTl framework has perhaps the widest user base with implementations of many formalisms, 
including stochastic process algebras (SPAs), stochastic automata networks (SANs) and generalised 
stochastic Petri nets (GSPNs). Mobius supports a distributed simulation environment and numerical 
solvers for models of up to tens of millions of states. 

PRISM f6] is a probabilistic model checker which supports low level formalisms such as DTMCs, 
CTMCs and Markov Decision Processes (MDPs) with an analysis engine based on Binary Decision Dia- 
grams (BDDs) and Multi-Terminal Binary Decision Diagrams (MTBDDs). PRISM can analyse models 
of up to lO'^ states, however this can depend heavily on the model being studied and on detailed consid- 
erations such as the exact variable ordering in the underlying MTBDD. 

Performance tools that support differential-equation based analysis have been primarily designed around 
stochastic process algebras such as stochastic Ti-calculus and PEPA. For TT-calculus SPiM lfT2l [T3]| has 
long been the standard tool for simulating stochastic n calculus models, but being a simulator it suffers 
from scalability issues for models with very large populations of components. A recent tool, JSPiM |[T4ll . 
allows for the chemical ground form subset of stochastic TT-calculus to be analysed via differential equa- 
tions. For the stochastic process algebra PEPA, the tools ipc fTS^.H^ and the Eclipse PEPA plug-in ITTI 
implement the so-called ^m/J translation [3] to produce sets of differential equations for PEPA. 

In the field of biological modelling, tools such as Dizzy ifTSll and SPiM have been used to capture first 
order approximations to system dynamics using a combination of stochastic simulation |[T9l and differ- 
ential equation approximation. A recent tool from GOl generates ODEs approximating higher moments 
of models using the mass-action kinetics and described in the Systems Biology mark-up language. How- 
ever for the reasons discussed, these techniques do not extend to computer and communication system 
modelling. 

The tool presented here provides a framework implementation for a variant of PEPA known as Grouped 
PEPA or GPEPA {2T\ . Using the established relationship between the underlying discrete Markov pro- 
cess and the fluid model [IJ, we can for the first time produce higher moments through differential 
equation analysis of massive models. The relationship with the underlying discrete model also allows us 
to identify areas of possible inaccuracy where the fluid model is known to be least accurate. Finally the 
GPEPA analyser (GPA) allows comparison by stochastic simulation, where this is feasible. 
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1.3 Grouped PEPA models 

Grouped PEPA or GPEPA [1] is a simple syntactic extension of PEPA. A brief summary of PEPA is 
given in Appendix |A] GPEPA is defined to provide a more elegant treatment of the ODE moment ap- 
proximation (and also allowing a more elegant implementation). Formally, the extension introduces a 
further level in the syntax of PEPA, the Grouped PEPA models, defined as: 

M::=Mt<iM \ Y{P\\ ■■■\\P} 

This defines a GPEPA model to be either a PEPA cooperation between two GPEPA models (over the set 
of actions L) or alternatively a labelled grouping of PEPA components, P, in parallel with each other. 
Y is the group label. The Grouped PEPA model is nothing more than the standard PEPA model with, 
additionally, a label to define the components involved in parallel grouping. These labels will be later 
used to define the level at which the ODE approximation to the system is made. 

The example from Section [TTT] can be represented by the following GPEPA definition: 

Processor^) = {acquire, ri). Processor \ Resource^ = {acquire, r2).Resource\ 

Processor\ = {task, q).Processoro Resource\ = {reset, s).Resourceo 

System = Processors{Processoro[m]} ^ ^ Kesources{ResourceQ[n\} 

The choice of labels (Processors and Resources) was arbitrary and will serve to identify each state of 
the system by counting the occurrences of Processor^ and Processori processes in the group Processors 
and occurrences of Resource^ and Resource \ in the group Resources. 

1.3.1 ODE analysis 

The states in PEPA models originally keep track of the state of each individual sequential component. 
This can lead to state space explosion, which makes the model not amenable to traditional analysis 
methods other than the computationally expensive stochastic simulation. The state space explosion is 
especially severe (with respect to the syntactical size of the model) in the case of models with groups 
consisting of many components acting in parallel. An established way to tackle this in the case of groups 
consisting of many identical components is by aggregating the state space by keeping track of counts of 
the individual components [3 |. In the context of Grouped PEPA models, it is sufficient to represent each 
state of the underlying CTMC by a numerical vector N(f) consisting of counts NG.p{t) for each possible 
pair of group label G and component P (as in |[T1). 

Allowing this to be a real-valued vector y{t) (with components accordingly named VG,p{t)), a system 
of ordinary differential equations, y{t) = f(v(f)) with initial conditions v(0) = N(0) can be intuitively 
derived from the corresponding PEPA description. These ODEs deterministically approximate the PEPA 
model's evolution over time. This has been formally shown in lUl (by comparison with the Chapman- 
Kolmogorov equations for the underlying CTMC) to be an approximation to the expectation of the indi- 
vidual group/component counts, i.e. y{t) f« E[N(f)]. 

The authors of fV\ further extended this method to derive ODE approximations to higher and joint mo- 
ments of the counts, such as variances, covariances and others. The number of differential equations 
generated by this method, when calculating all the moments of order up to p, is proportional to n^ where 
n is the number of component derivatives. 
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The method generahses the vectors N and v to include components for values (integers and reals respec- 
tively) of each possible moment of the group/component pair counting processes. For example N(f) can 
include the squared count of component P within a group G, and \{t) an approximation to its 

expectation, i.e. V(p^c)2(?) ^ E[A^jp^^(f)]. 

A general advantage of this approximation to the moments (including the mean) is that the resulting sys- 
tem of ODEs can be numerically integrated by standard methods. This usually requires less computation 
than running sufficiently many replications of the respective stochastic simulation over the same model. 



1.3.2 Two-stage Client/Server Example 

We will look at a more complicated example that will be used in the following investigations. Con- 
sider a client/server model, where communication occurs in two stages. Clients first request a server to 
communicate with. After admitting a request from a server, a client moves to a waiting state and the 
corresponding server to a ready state. Then, any ready server can serve a waiting client, after which the 
client can perform a required task and the server returns to the idle state. Additionally, an idle server can 
break, requiring a reset. This can be represented by the following PEPA components: 

Client = {request, r^eq) -Client ^waiting 
Client .waiting = (data , r^ata ) • Client Jhink 
Client Jhink = {think, Vthink) -Client 

Server = {request, rreq) -Server _get + {break, rtreak) -Server Jbroken 
Server _get = {data, rdata) -Server 
Server Jjroken = {reset, rreset) -Server 

composed into the GPEPA system equation defined by CS{c, s) over the set of actions L = {request, data} : 

CS{c,s) = C^tn\s{Client[c]] \^StryQrs{Server[s]] 

The aggregation described above represents each derivative state of the model at each time f as a vector 
N(f), i.e. the underlying CTMC can be treated as a vector valued stochastic process, with components 
^ciients.p(0' where P can be Client, Client ^waiting or Client Jhink and A'servers.gCO where Q can be 
Server, Server_get or Server Jbroken. The ODE approximation generates a system of ODEs, y{t) = 
f(v(f)) with a real-valued vector solution y{t) with the same components as the vector N(f) and initial 

conditions Vciients,aienr(0) = A^CIieiits,C«e«r(0) = C, Vservers,Sen.er(0) = A'^Servers.Sen'erCO) = S and Vg,p(0) = 



Otherwise. See Section 1.5 for the exact formulation of f(v(f)). 



1.4 Nature of the approximations 

The nature of the approximation of the system of ODEs to the CTMC evolution boils down to the ap- 
proximation of the rate at which two components evolve in cooperation (as discussed originally in ifTl). 

This is trivial in the unsurprising case of purely concurrent models, where no cooperation takes part. For 
these, the approximation is exact and can directly replace moment calculation via stochastic simulation. 
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In general, the rate of cooperation can be of the form r(N(f)) • min(f(N(f)),g(N(f))) where r(-) is a 
rational function of piecewise linear functions and f(-) and g(-) are piecewise linear functions all of the 
components of N(?). The functions f(N(f)) and g(N(f)) are a calculation of the total rates of cooperating 
actions being enabled in cooperating component groups. 

For a certain class of models called the split-free models, introduced in [jj, it can be shown that the 
function r is constant and f , g are piecewise linear. It then turns out that the nature of the deterministic 
approximation of N(?) by v(f) depends on the approximation 



E[min(f(N(0),g(N(0)] ~ min(f(E[N(0]),g(IE[N(0])) (1-1) 

where the right hand side is approximated by min(f(v(f)),g(v(f))). Note that the functions f and g 
may also include further minimum terms themselves and thus induce multiple further applications of the 
approximation not shown explicitly above. It is argued in [ 1] that the error of this approximation is at its 
highest around so-called switch points. These occur when the total rate of cooperating action between 
component groups becomes equal and hence f(N(?)) = g(N(?)) causing the min function to switch. 

One of the main contributions of this paper is to verify and further investigate this claim empirically. 
We will be able to use our GPEPA tool to identify the switch points and to display the error in the ODE 
approximation around the switch points of a given GPEPA model. We will focus on first order switch 
points: those coming from a min term involving no higher orders than means of component counts. 

The switch points are also defined for general GPEPA models. However, for those, the nature of the ODE 
moment approximation also depends on the approximation of terms, E[r(N(?))] for a rational function 
r(-). Our GPEPA tool allows analysis of these models, by using the approximation r(E[N(f)]). However, 
in our first investigation we concentrate on the split-free models. We will be using the client/server model 
CS{c,s) from the above section, as it can be easily shown to be split-free. 



1.5 Theoretical considerations 

In this section, we aim to provide justification for the convergence of the suitably-scaled mean and 
variance ODE approximations as the component population sizes are increased. 

In order to do this, we will more formally set up some notation to allow us to consider general GPEPA 
models and their associated systems of ODEs in a compact manner. For a general GPEPA model, assume 
that we have a vector- valued stochastic process N(f), defined on M+ taking values in C M+, for some 
G Z+. In line with Section 1.3. 1[ each component of this process counts the number of a particular 



derivative state currently active in a parallel group of the model, of which there are N derivative states in 
total, across all parallel groups. 

Analogously, we define the vector-valued deterministic function v(f), also defined on and taking 



values in to be the ODE approximation to the expectation of this model. We assume that it is 
defined uniquely by the system of differential equations \{t) = f(v(f)) and the initial condition v(0) = 
N(0). As discussed above, the exact form of the deterministic function f : M.^ — M'^ is given in fill 
and corresponds component-wise to the rate at which each derivative state is incremented, minus that 
at which it is decremented, in a given state of the model. It helps now to consider the system of ODEs 



explicitly for an example model. Indeed, in the case of the client/server example of Section 1.3.2 we 
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have a total of six derivative states, giving N = 6, 



v(0 



{Nc{t),Nc„{t),Nc,{t),Ns{t),Ns^{t),Ns,M'' 
M), vcAt), vc,(0, V5(0, vs,(0, vsM'^ 



and: 



-mm{vs{t), Vc{t))rreq + Vc,{t)rthmk 
-min(vc„(0, vsg{t))rdata + v[nri{vs{t), vc{t))rreq 
-vc,{t) rthink + min ( vq, {t),vs^{t))rdata 

mm{vs{t), Vc{t))rreq - Vs{t)rbreak + min(vc„„ [t], Vs^ {t))rdata + VS,, {t)rre.set 

-min(vc„,(0) vs^{t))rdata + v[m^{vs{t), Vc{t))rreq 

-VSb {t ) Preset + V5 (0 rbreak 



f(v(0) 



v 



where we have used an appropriate shorthand notation for the components of N(f) and \{t). In order 
to make the desired theoretical considerations in this section, it is necessary to consider a structurally 
equivalent sequence of GPEPA models which have the same structure of parallel component groups and 
differ only in terms of the size of the component populations within these groups. Furthermore, we 
require that they all have the same initial proportion of each component type in each case. For such 
a sequence of GPEPA models, let {N'(f)}Jlj be their associated stochastic counting processes in the 

notation above, and for each /, write Si := N[{t) -\ +NI,{t) for the total component population of 

model irjSo our requirement of constant initial component type proportions is stated formally as: 



In the case of (G)PEPA, it is relatively straightforward to see that for any x G M^, f(^x) = ^f(x) for 
all k G R+. Furthermore, since the GPEPA models in our sequence differ only in terms of their initial 
component counts, it is easy to see that the function f(-) is the same for any /. These two facts together 
mean that we need only define the fluid approximation to N'(f), say, \'{t) for a particular value of 
/, and the fluid approximation for any other / can be defined in terms of it. Indeed, for any /, j > 0, if 
v'(f) = f(v'(0) and v-''(f) = f(v-''(f)) with initial conditions, v'(0) = N'(0) and v-''(0) = N-^^O), respectively. 



Thus for the rest of this section, we consider the quantity \{t) := \'{t)/Si, for all f > 0, which we have 
just seen is independent of /. 

1.5.1 Convergence of mean approximation 

It is shown in lITOll . based on a result by Kurtz |[22l . that, in the above notation for any sequence of 
structurally equivalent GPEPA models, the following holds: 
Theorem 1.1 If Sj — ^ oo as / — > oo then: 



N'{0)_Nl{0) 



for all /, 7 > and k £ {I, ...,n} 



S,- Sj 



we have: 



v'(0) = v^'(O) X — and y'{t) = yj{t) x — for all f > 





as i — > oo, uniformly over bounded intervals t G [0, T). 



^This does not depend on t because the PEPA operational semantics preserve component populations. 
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This demonstrates that, in the limit of component population size, the mean number of components 
at time t does indeed tend towards the ODE solution. An example of this convergence is shown in 
Figure [6j a). Section [3] 



1.5.2 Convergence of variance approximation 

In this section, we provide theoretical support for our hypothesis that in the limit of large populations and 
under an appropriate scaling, the variance of component counts converges to the approximation given by 
the ODEs for spht-free GPEPA modelsj^ 

We will start by approximating the GPEPA model's underlying CTMC, N'(f), as the sum of a determin- 
istic process, \{t), given by the first order ODEs, and a Gaussian process, E(?) defined below. From this 
description, we will derive a set of ODEs describing the evolution of the covariances of the Gaussian pro- 



cess. These can be formally shown (Theorem 1 .2 1 to agree with the covariances of the CTMC in the limit 



of a sequence of GPEPA models of increasing size. We will argue that the second moment ODEs from 



the CTMC (Section 1.3.1 1 capture the dynamics of the system more accurately than the ODEs generated 



from the Gaussian process. This provides the basis for our hypothesis that the variances of the CTMC do 



indeed converge to the solution to the second moment ODEs from Section 1.3.1 as the population size 
increases. 

The decomposition described above gives the following approximation to the underlying CTMC of a 
GPEPA model: 



N'(0~S,-v(0 + V'^/E(0 (1.2) 

where E(f) is the Gaussian process mentioned above. In order to proceed with defining E(f), it is neces- 
sary to enumerate explicitly the transitions in the aggregated state space of the GPEPA model. Assume 
that there are K such transitions and, following |[23l . let {1^ € Z'^jf^j be the sequence of jump vectors, 
specifying that if the ^th such transition occurs at time t, N(f) = N(f— ) +1*^, where t— is the instant 
immediately preceding t. Then define rate functions, {f^ : — )• M+}f^j, specifying the aggregate rate 
of each transition. In the case of the client/server example, we have K = 5 (for the 5 possible activities) 
and: 

1^ = (-1, 1,0, -1, 1,0) f{x) = mm{xi,X4)rreq 

l2 = (0, -1, 1, 1, -1,0) /2(x) =mm{x2, X5)rdata 

1^ = (1,0, -1,0,0,0) f{x)=X3rMnk 

1^ = (0,0, 0,-1, 0,1) /^(x)=X4W 

1^ = (0,0, 0,1,0,-1) f{x)=X(,rre,et 

The stochastic process E(f) is now given in the following theorem ||23l . which applies to any sequence 
of structurally equivalent split-free GPEPA models with underlying CTMCs {N,(f)}^j, corresponding 
model sizes and rescaled ODE approximation \{t). 

Theorem 1.2 below defines the Gaussian process E{t) in terms of time-scaled Wiener processes. 
Readers unfamiliar with Wiener processes and weak convergence of stochastic processes can consult for 
instance Rogers and Williams Il24l or Klebaner |[25l . The set T defined in the theorem can be given more 
intuitively as the set of all times at which two arguments of a minimum function in the right-hand side, 
f(v(f)), of the first moment ODEs are equal. 



Similar considerations should also be possible for non split-free models, but we do not consider these here for brevity. 
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Theorem 1.2 Let T > and let t be the subset of {t £ [0, T)} for which f(-) is not totally differentiable 
at the point \{t). We require that T has Lebesgue measure zero. Then on all of [0, T) \ T, f(-) has a 
well-defined Jacobian at the point \{t), say D({\{t)). Extend this to all points {\{t) : ? G [0, T)}, say by 
defining it to be the matrix of zeros at times in T. 

Then if S; — )■ oo as / — )• oo, — ^/Si\{t) E(f) aiso as / — )• oo, where: 

E{t) := rDf(v(s)) •E(.)d.+ X W' ( ffiy{s))ds) 

Jo J 

and 

J is a sequence of K mutually independent standard Wiener processes (aka Brownian 
motions). The convergence (^) is weak convergence in D^n[0, T), the space of M.^ -valued cadla^' 
functions. 



This is essentially a generalisation of a result of Kurtz 1*261 to cope with the case of non-smooth rate 
functions as introduced in the case of PEPA models, or other formalisms for modelling synchronisation 
in computer systems, by the minimum functions. It can be shown that the process E(f) is the unique 
solution of the following (Ito) stochastic differential equation (SDE): 

dE{t)=n{E{t),t)dt + a{t)dW{t) 

where ^(x, ?) : x M+ ^ and o{t) : M+ M.^""^ are defined by: 



M(x,0:=Z)f(v(0)-x, 



a(0 := (//x^/i(v(0) 



and W{t) := {W^{t), . . . ,W'^{t))'^ is a ^-dimensional standard Wiener process. This representation 
allows us to apply the machinery of Ito's Lemma, see for example ll24l l27l . to derive the following 
system of ODEs, whose unique solution is exactly the covariance matrix of E(f): 



-Cov[E(0, E(0] = Cov[E(0, E(0] (Df(v(0))'^ + Df(v(0) Cov[E(0, E(0]"^+ £ /(v(O)l^ (1^)'^ 

keK 

(1.3) 

If we apply this to the client/server example for the specific covariance Es{t), corresponding to the Server 
component of N(f), Ns{t), we have: 

^Cov[Es{t),Esit)]= -2r,,^(l{,^(,)<,^(,)}Cov[£5(0,'Es(0] + l{v,w>vcW}Cov[£5(0,'^c(0]) 

+ 2rdata (l{v,Jr)<vc,^,(f)}Cov[£5(0> Es^{t)] + ^{vs,{t)>vcJt)}Cov[Es{t), Ec„{t)] 
-2rbreakCo\[Es{t),Es{t)]+2rresetCo\[Es{t),Es^{t)] 

+ rreqrmn{vc{t), + ''datomin(vc„,(0, Vs^{t))+rhreakVs{t) + rresetVs,{^) 

(1.4) 

Note that Theorem [L2] suggests the approximation: 

Coy[Ns{t),Ns{t)\ « Coy[Svs{t) + VSEs{t), Svs{t) + VSEs{t)] = SCoy[Esit), Es{t)] 



^Continue a droite, limitee a gauche, that is, right continuous with left hmits. 
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where S is taken as the total population of clients and servers. Indeed, in general. Theorem 1 1 .21 implies 
(assuming its hypothesis), as / — )• oo, if Si — oo, that: 

icov[N'(0, N'(0] ^ Cov[E(0, E(0] 

So the system of ODEs given in Equation ( |1.3| ) yields an approximation to the covariance matrix of the 
underlying CTMC of a general GPEPA model. Furthermore, as we have just seen, this approximation is 



guaranteed by Theorem 1.2 to converge in the limit of large populations. 

However, the covariance ODE approximation implemented by the GPA tool consists of integrating a 
slightly different system of ODEs. Our intention now is to show that they are very similar to those of 



Equation ( 1.3 1 and in fact can intuitively be regarded as a better approximation to the actual covariance 
matrix of the CTMC. This is the basis of our conjecture that a similar convergence result also holds, and 
furthermore, that the rate of convergence may well be faster for the ODEs implemented in the GPA tool. 

In more detail for our example, applying the techniques of UJ for the specific term Cov[Ns{t), Ns{t)], we 
obtain the exact differential equation: 

^Cov[Ns{t),Ns{t)] = ^E[Nl{t)]-2E[Ns{t)]^E[Ns{t)] 



- 2r,,, (^E[min(A^c(0 A^5(0, A^i(O)] - IE[min(Afc(0, ^5(0)]E[A^s(0] 
+ 2r,ata (E[imn{Nc,,{t)Ns{t),Ns^Nsm - E[min(A^c„,(0, A^s,(0)]E[A^5(0] 
-2rbreak(nNl{t)]-E^[Ns{t)] 



+ 2 w l^nNs, it)Nsit)] - E[Ns, (t)] E[Nsit)] 
+ rre,nmmiNcit),Nsm+rdatanrnm{NQXt),Ns,m 

+ rbreaknNs{t)]+r,-ese,nNs„{t)] (1-5) 

Since the corresponding system of ODEs cannot be solved analytically or numerically due to the pres- 
ence of expectations of non-linear functions, the approximation E[min(X, Y)] ^ min(E[X], E[F]) can be 



applied repeatedly as in HI and Section 1.4 to deduce approximate ODEs for Cov[Ns{t), Ns{t)] and the 
other covariances. This system of ODEs can then be solved numerically as implemented in general by 
the GPA tool. 

There are two kinds of approximations being applied here, those which we might call first-order, such 
as K[mm{Nc{t) , Ns{t))] « min(E[A'c(?)]7 ^i^sit)])^ and those which we might call second-order, such 
as E[mm{Nc{t)Ns{t),Nj{t))] ^ mm{E[Nc{t)Ns{t)], E[A^|(0]). The key point to note now is that if we 
keep the first-order approximations the same but replace the second-order approximations by ones of the 
form: 

E[min(A^c(0A^5(0,A^I(0)]~l E[Nl{t)] 

then we recover the system of ODEs of Equation ( |1.3| ). This is a reasonable approximation, but we 
switch second-order minimum terms by making only first-order comparisons. It is thus intuitively clear 
that it is likely to be a worse approximation than the second moment ODEs derived from the CTMC. 
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On this basis, we might reasonably expect that the covariance ODE approximations, derived from the 
CTMC and implemented by the GPA tool, should converge to the actual covariances when scaled by the 
total component population size. An example of this can be seen in Figure[6jb), Section[3] 



2 GPA: The GPEPA Analyser 



In this section we introduce a new tool for analysing Grouped PEPA models. The Grouped PEPA anal- 
yser (GPA) uses the results from |[T1, briefly described above, as a basis for producing deterministic 
approximations of transient behaviour of syntactically specified Grouped PEPA models. It can also use 
these to provide passage-time approximations from lITOll . 

GPA in addition implements stochastic simulation of the models to allow investigations into the nature 
of the approximation and the convergence results from Section 1.5 In the next section, we look at some 
specific examples that empirically demonstrate this. We first give an overview of GPA, giving its syntax 
and the commands it provides. 



2.1 Overview 

The syntax of models specified in GPA is very close to the formal language of Grouped PEPA. Each 
GPA model starts with definitions of parameters used as rates in component definitions and parameters 
used as initial counts in the grouped model definitions. The definitions of individual PEPA components 
then follow. Finally, a single Grouped PEPA model is given, using the defined components. 

Groups are specified by labels with components enclosed in braces { }. Multiple components of the same 
type are given by [n] written after the component identifier, where n is a previously-defined parameter. 
The cooperation operator is represented by the cooperation set enclosed in angled brackets < >. 

The GPEPA model corresponding to the example from Section [TTT] can be represented in GPA as: 



r 1=2.0; q=14.0; m=50.0; r2=14.0; s=2.0; n=20.0; 

ProcessorO = (acquire, rl) .Processorl; ResourceO = (acquire, r2) .Resourcel; 
Processorl = (task, q) .ProcessorO; Resourcel = (reset, s) .ResourceO; 

Processors{ProcessorO [m] }<acquire>Resources{ResourceO [n] } 



Appendix |C] contains the grammar describing the full syntax of GPA. 
2.2 ODE Analysis and Comparison with Simulation 

On each GPA model, several analyses can be performed. Each analysis is specified with required param- 
eters (for example, the time range we are interested in for the transient behaviour of the model), after the 
analysis name. Following this is a list of commands that allow visualisation and further manipulation of 
the resulting data points. 

The ODE analysis provides analysis of the given model by a deterministic approximation via a set of 



ordinary differential equations as described in Section 1.3.1 The enclosed commands determine the 
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maximum order that has to be considered in calculations. For example, if the commands involve only 
plots of means and variances of component counts, differential equations for approximation of all the 
possible (joint) moments of order 1 and 2 are implicitly generated. These ODEs are then numerically 
integrated (using a fourth-order Runge-Kutta solver). The parameter stopTime specifies the time over 
which the numerical solution is given. Parameter stepSize determines the fixed time interval between 
each successive pair of data points that are taken. The parameter density specifies the accuracy by 
telling the solver how many sub-intervals should each time interval between the data points be divided 
into. The following GPA code performs an ODE analysis on the above GPA model, displaying: mean 
counts of Processoro components in group Processors and Resource^ components in group Resources 
and model switch points: 



odes (stopTime = 3.0, stepSize = 0.001,d.ensity=10)-[ 

plot (E [Processors : ProcessorO] , E [Resources : ResourceO] ) ; 
plotSwitchpoints (1) ; ... 

} 



The Simulation analysis provides analysis of models by stochastic simulation of their underlying con- 
tinuous time Markov chain. GPA generates a representation of the CTMC and uses the Gillespie algo- 
rithm ||T9 1 to simulate the CTMC at given time intervals until the simulation stop time is reached. 

The stochastic simulation is repeated several times (given by the parameter replications) and then 
averaged to provide the final data set. On the above model, the following code performs a simulation 
analysis to extract the same mean component counts as before: 



simulation(stopTime = 3.0, stepSize = 0.001,replications=1000){ 
plot (E [Processors : ProcessorO] , E [Resources : ResourceO] ) ; . . 

} 



To compare results from stochastic simulation and ODE analysis, the Comparison analysis provides a 
way of calculating the difference between the two analyses. It takes a simulation and ODE analysis as 
parameters (which are required to have the same stop time and step size) and calculates the difference 
between their resulting data points. All the enclosed commands then use this difference in the same way 
as data points for the other analyses. For example, a comparison of the above analyses: 



comparison ( 

odes (stopTime=3.0,stepSize=0. 001, density=1000){. . .}, 
simulatioii(stopTime=3 . , stepSize=0 . 001 , replications=1000) {...}){ 
plot (E [Processors : ProcessorO] ,E [Resources : ResourceO] ) ; 

} 



2.3 Commands and Functionality 

GPA can plot results from these analyses itself or output raw data by means of an optional file redirection. 

The plot command provides direct plots of different arithmetic expressions involving (higher order and 
joint) moments of component counts and numerical constants. Each group-component pair is identified 
by the syntax group : component. A moment is an expectation, with syntax E [] , of products of these. 
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for example E [Gl : CI "2 G2 : C2] . Several shorthands are provided, such as the variance represented by 
Var [G : C] (which stands for E [G : C"2] -E [G : C] ~2), and Gov [Gl : CI , G2 : C2] which represents the co- 
variance of two component count variables and the «th central moment represented by Central [G : C , n] . 
The following commands are examples of the above: 



plot(Var [Processors :Processorl] , Var [Resources : Resourcel] ) ; 

plot (Central [Processors : Processor 1 , 4] ) ; 

plot (E [Resources : ResourceO] +E [Resources : Resourcel] ) ; 



The plotSwitchpoints command inside an ODE analysis visualises distance from switch points in the 
ODE approximation. It obtains a set of all the occurrences of the min function (containing only moments 
of order upto a given integer argument) in the moment ODEs. For each data point from the ODE analysis, 
the command plots the difference between the two arguments of each of the min occurrences. The switch 
points then correspond to the zero points of this plot. 

One of the reasons the plot command provides arithmetic expressions of the moments is to give GPA 
the flexibility to obtain approximations of passage times. As described in |10], the passage time ap- 
proximations and the corresponding bounds can be expressed as certain functions of the (higher order) 
moments. We will explore this functionality in a later paper, but suffice to say that the accuracy of the 
passage-time approximations depends critically on the ODE approximation of first and higher moments. 

Some example plots showing error comparisons in variances for the processor/resource model can be 
found in Appendix |B] For further details on GPA, including random model generation for testing pur- 
poses, we refer the interested reader to the freely available source code obtainable from [,28,1 . 



3 Numerical Investigation: Two-stage Client/Server 

In this section, we investigate empirically the nature of the ODE approximations as mentioned in Sec- 



tion 



1 .4 and the convergence results from Section 1.5 More specifically, we will check that the simulation 



variances converge to the ODE approximations as the initial component populations scale up. We also 



look at higher moments to suggest a possible result analogous to the one mentioned in Section 1.5.2 



We will look at the eiTor of the ODE approximation of the moments and quantitative properties of the 



convergence, in the context of the switch points defined in Section 1.4 

In order to investigate this, we use two models - one which does not stay near a switch point in any 
large time interval (we will informally refer to it as to the occasionally switching model) and one which 
steadily stays near a switch point for a longer period of time (we will informally refer to it as to the 
persistently switching model). 



In both cases we will use the two-stage client/server model from Section 1.3.2 with two sets of parame 



ters and taking the client population, c = 100, and the server population, s = 50: 

Model A: r^^^ = 2.0 rhreak = 0.l r,ft,v,i = 0.20 r^ato = 1-0 r,.eset = 2.0 

Model B : r^^ = 2.0 rbreak = 0.3 rMnk = 0.35 r^ata = 2.0 rreset = 0.05 

The switch point behaviour for these examples can be seen on Figure[4] There are two possible sources of 
switch points in the client/server model, each corresponding to an instance of cooperation in the model. 
One, the term min(vc/;^-„r(0 • recj 1 ^Server (0 • rreq), comcs from the cooperation when a client estabhshes 
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connection with a server. Another, mm{v aient,,,,,,,,,^ (0 • rdata,vserverg„ (0 • tdata), comes from the cooperation 
when a client retrieves data from an available server. 

For the min term involving Servergf.t and Client ^^,aiting, the model hits infinitely many switch points. 
These do not influence the error of the approximation since the two corresponding counting processes 
are stochastically identical, so there is no error in the corresponding expectation approximation. For the 
min term involving the Client and Server components, in the first, simpler case, the model hits one switch 
point at time f ~ 2.1, but does not stay around any switch point for a period of time. The second model 
hits two switch points when f ~ 2.8 and ? ~ 4.8 and stays close to a switch point in the time interval 
between them. 




(a) Switch point plot for Model A (b) Switch point plot for Model B 

Figure 4: First order switch points for tlie two-stage client/server models 



In the following sections, we will look at the errors in moment approximation for both of the above 
models and compare them in the context of these switch points. We will look at the expectation, variance, 
skewness and kurtosis (the standardised third and fourth central moments) of the component counts for 
each model and its respective versions with initial populations scaled by a factor of n = 1,4, 16 and 64. 
We investigate how the scale influences the error of the ODE approximation. For each scale, we plot 
the difference between the moments from simulation and their approximation, specially near the times 
where a switch point is encountered. It is worth noting that these times are the same across all n, since 



the ODEs for the means are scale invariant, as was mentioned in Section 1.5 



In line with the considerations from Section |1.5[ we normalise the error of expectation and variance 
approximations dividing by the total component population. We plot the errors of the skewness and kur- 
tosis approximations (standardised third and fourth central moments, respectively) without any additional 
normalisation. 



3.1 Model A: Occasionally switching model 

First, we look at the simpler case. Model A, in which the client/server model does not stay near a switch 
point for a longer period of time. We start by comparing side by side, for the first three central moments, 
the simulation average with its ODE approximation in Figure [5] 

For kurtosis (not shown), similar to the skewness, the approximation is quite accurate when the model is 
not near a switch point, but errors are visible otherwise. 
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(e) ODE approximation: skewness (f) Simulation: skewness 

Figure 5: Moments from simulation and approximation for the two-stage client/server model, Model A. 



For all four moments the ODE approximation is qualitatively close to the simulation average. However, 
there are visible quantitative differences starting at the variance and getting worse with the skewness and 
kurtosis. This is especially noticeable for the moments involving the number of Client components and 
in the time interval around the first order switch point from Figure |4] (shown by the vertical line on the 
plots in Figure |5]l. We plot the actual error of the approximation (using the GPA comparison analysis) 
for these moments, for each of the scaled versions of the model, in Figure [6] 
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3.2 Model B: Persistently switching model 



Figure [7] shows the mean and variance for the case where the model stays near a switch point for a longer 
interval of time in Model B. The mean seems quite accurate, but the error in variance approximation is 
high around the interval where the model stays near a switch point. 

Figure [8] looks at the error more closely and plots the difference between the moments from simulation 
and their ODE approximation for the Client component, for different scales of Model B. It can be seen 
that the normalised error in Figure [8] is much higher than in the case of Model A. We believe this is 
caused by the fact that the model is closer to the switch points. However, in both cases we can confirm 
that the error for the mean and variance seems to be going to zero in the scale limit, but concentrated 
most around a switch point. The same seems to be the case for the skewness. For kurtosis, it seems that 
the error does not necessarily get smaller in value, but the interval where it appears seems to get smaller 
with increased scale. 



Section 1.5 shows that the marginal distributions of component counts become normal as the populations 
are increased. The normal distribution has a constant skewness and kurtosis (0 and 3 respectively). We 
can check the plots for the simulations to verify this. This seems to be the case for the ODE approxima- 
tion, except perhaps for the points concentrated around the switch point times. 
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(c) ODE approximation: variance 
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(d) Simulation: variance 



Figure 7: Moments from simulation and approximation for the two-stage client/server model, Model B. The 
ODE approximation is quite accurate for the mean, but large differences are visible for the variance. 



4 Conclusion and future work 

We introduced a tool, GPA Il28l . that for the first time makes available the ODE approximations of higher 
moments to a wide range of models described in an extension of the PEPA stochastic process algebra. 
We used the tool to carry out investigations into the nature of the ODE approximation, that will help the 
modellers to detect errors without running computationally more expensive stochastic simulations. We 
theoretically justified that the variance of the component counts converges to the ODE approximation as 
the initial populations are scaled and with the help of our tool verified this for a simple example. 

We observed, for a model where the resulting differential equations are piecewise linear, that the error 
is influenced by how closely the model stays near the switch points during its time evolution. If the 
model only crosses switch points at certain points of time and does not stay near any during the rest 
of the time, then the error is concentrated tightly around those switch point events. We also saw that, 
with increasing scale of initial component populations, the error in the ODE solution becomes even more 
tightly concentrated around the switch point. If the model stays near switch points for longer periods 
of time, the resulting error is much more severe and decreases more slowly with increased scale of the 
initial populations. 

These observations help us to assess the validity of the ODE approximation without actually running the 
simulations. For a given model, we can use GPA to visualise the switch point behaviour of the model 
and use the intuition gained from our investigations to say whether the approximation is accurate. 
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Moreover, the presented results and observations are not just specific to the PEPA stochastic process alge- 
bra. The min functions with the concept of switch points appear in situations when there is a competition 
for multiple resources, for example multi-server queues or many-server semantics stochastic Petri nets. 
Therefore we believe that the gained insight is relevant to a wider area within performance modelling of 
computer systems. 

In future, we plan to develop methods that would be able to quantify the error (e.g. in terms of bounds 
obtained from the distance from a switch point), thus making GPA able to warn the modeller of the 
potential magnitude of any error. We also plan to investigate the sensitivity of the switch point behaviour 
to changes in the rates and initial population parameters to allow deliberate generation of accurate ODE 
approximation. 
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A PEPA Summary 

We briefly introduce PEPA f29\ , a simple stochastic process algebra with sufficient expressiveness to 
model a wide variety of systems. As in various other process algebras, systems in PEPA are represented 
as compositions of components which undertake actions thus evolving into further components. Each 
action has a rate, interpreted as a parameter of an exponential random variable that governs the delay 
associated with the action. This means that the stochastic behaviour of the model can be described as a 
CTMC. 

The basic building blocks of PEPA are sequential components, described by the syntax 

S::= {a,r).S \ S + S \ Cs 

where Cs stands for a constant that defines a sequential component. Intuitively, the component {a,r).S 
can undertake the action a with rate r and evolve into the component 5. The component P + Q has a 
choice to undertake actions that both P and Q can undertake. 

The sequential components can be composed in parallel to form model components, described by the 
syntax 

P : =p\X]p \ P/L \ C 

where C is a constant defining any PEPA component. The component P IXl g enables cooperation be- 
tween P and Q by only allowing P and Q to evolve together when undertaking an action from the co- 
operation set L. The cooperation influences the rate of the common evolution. PEPA assumes bounded 
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capacity: a component cannot be made to perform an activity faster by cooperation, so the resulting rate 
is the minimum of the cooperating rates. This is discussed in more detail in [29 1. For actions not in L, P 
and Q can evolve independently with no influence on the rates. If L is the empty set, we use the shorthand 
P II 2 meaning that P and Q are purely concurrent and don't synchronise. We also use the shorthand P[n] 
for n purely concurrent copies of P. 

The syntax P/L stands for action hiding. For simplicity, we will not consider this feature in the presented 
work, however, it should be straightforward to include using the results from HI. 

The full details of the operational semantics of PEPA can be found in ||29ll . 

B Processor/Resource Variance Plots 




(a) Normalised difference in variance (b) Normalised difference in variance for the scaled version 

Figure 9: Difference in variance between the simulation and its ODE approximation for the proces- 
sor/resource model, produced by the comparison analysis. 

Figures [T] [2] and [3] that we used in the initial example in Section |1.1| were all produced by the respec- 
tive GPA commands enclosed in simulation and ODE analyses. We could examine the error of the 
approximation more closely by using the comparison analyses. Figure |9] compares the error in variance 
approximation for the processors resource model with m = 50 and n = 20 with a version with populations 
scaled by 100, i.e. with m = 5000 and 2000. In order to demonstrate the results from Section [T3| the 
difference is scaled by the population size m + n. We can see that under this scaling, the relative error 
does decrease. 
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C GPA syntax definition 

Models 

System := ParameterDefinition* ComponentDefinition* ModelDefinition Analysis* 

ParameterDefinition := parameterld = realnumber; 

ComponentDefinition := componentID = Component 

Component := Component<ActionLisf >Component 

I Summation \ componentid \ (Component) 
Summation := Prefix{+ Prefix)* 

Prefix := {actionid, parameterld ) .{{Summation) \ stop | componentid |); 

ModelDefinition : = {ModelDefinition <ActionList^ >ModelDefinition) 
I groupLabel{ComponentsParallel} 
ComponentsParallel := Component{ \Component)* 

Component := componentid {[parameterld])^ 
ActionList := actionid {, actionid)* 

Analyses 

Analysis := ODEs \ Simulation \ Comparison 

ODEs := odes(stopTinie =realnumber, stepSize =realnumber, density =integer) 
{Command*} 

Simulation := simulation(stopTiine =realnumber, stepSize =realnumber,zeplica.tioiis =integer) 
{Command*} 

Comparison := coiiipa.rsioii{ODEs,Simulation){Command*} 
Commands 

Command := CommandN oFile{->" filename")^' ; 
CommandNoFile := plot{MomentExpressions) \ plotSwitchpoints(mteger) 
MomentExpressions := MomentExpression{, MomentExpression)* 
MomentExpression := MomentExpression{+ | — | * | / | ")MomentExpression 
I E[Moment{+Moment)*] \ Ma.-r[GCPair{+GCPair)*] 
I (Standardised) 'Central [GCPi3/r,mteger] 
I realnumber \ parameterld \ {MomentExpression) 
Moment := {GCPair{^integer)'')^ 
GCPair := groupLabehComponent 



